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SUMMARY 


“Fd2"  is  a  software  package  developed  at  Tetedyne  Geotech  Alexandria  Laboratories  (TGAL)  dur¬ 
ing  the  past  several  years  for  genersrtfog  synthetic  seismograms  and  displaying  the  wavefieids.  This 
package  consists  of  primarily  a  2-dlmort8lonal  2nd-order  expficil  inear  flnite-differorx»  (LFD)  code.  LFD 
method  has  the  advantage  that  the  solution  contains  all  conversions  and  an  orders  of  muttiple  scattering. 
It  permits  examinations  of  fairly  general  models  with  arbitrary  complex  variations  In  material  properties 
and  free-surface  geometry.  Furthermore,  it  does  not  require  many  assumptions  commonly  Invoked  bi 
other  theoretical  approaches.  The  basic  Imitations  to  the  LFD  method  or  the  finite-element  method  are 
the  computational  cost  and  memory  requirements.  These  constrain  the  size  of  the  grid  and  the  number 
of  time  steps  that  can  be  calculated  over  a  reasorwfole  time  frame.  Our  LFD  code  has  a  distinguishable 
feature  in  that  it  allows  the  inclusion  of  topogrj^lcal  free  surface.  This  is  particularly  useful  in  modeling 
nuclear  explosions  buried  in  mountains. 

In  this  topical  report,  sample  scripts  are  presented  to  illustrate  the  usage  of  'yd2"and  several  sup¬ 
porting  routines  for  plotting  out  the  synthetics,  generating  2-dimensional  media,  as  weU  as  the  graphic 
visualization  of  wavefieids.  The  algorithms  for  handing  the  boundary  conditions  of  polygonal  topography 
are  reviewed  in  detail.  Thus  this  topical  report  serves  as  both  a  programmer's  guide  and  the  user's 
manual. 
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1.  INTRODUCTION 


Wave  propagation  problems  in  seismology  involve  the  solution  of  a  set  of  differerttial  equations  in  a 
medium  in  which  the  material  properties  vary  in  spa<^,  i.e. .  the  earth.  The  use  of  numerical  simulations 
is  a  straightforward  means  for  studying  this  kind  of  problem  especially  when  laterally-varying  velocky 
stnjcture  or  complex  topographic  relief  is  encountered.  Methods  such  as  Gaussian  beam  technique  and 
ray  theoretical  schemes  are  restricted  to  cases  where  variattons  of  the  medium  are  much  larger  than  the 
seismic  wavelength.  The  Kirchhoff-Helmholz  integration  method  is  useful  for  media  with  sharp  inter¬ 
faces,  but  it  does  not  include  the  multipie  scattering  along  the  irteriaces.  and  it  is  not  appropriate  for 
reflections  from  velocity  gradients  similar  in  extent  to  the  seismic  wavelenglh.  Perturbation  methods  are 
applicable  only  for  weak  scattering  problems.  Among  all  numerical  approaches,  finite-difference  (FO)  and 
finite-element  (FE)  methods  are  not  restricted  to  velocity  variations  of  a  particular  size  with  respect  to 
wavelength.  FO  and  FE  can  generate  synthetic  seismograms  for  very  compBcated  media  in  cases  of 
weak/strong  or  multiple  scattering. 

FD  method  solves  either  the  wave  equations  or  the  elastodynamic  equations  by  replacing  the  par¬ 
tial  derivatives  in  space  and  time  by  their  finite-difference  approximatiorts.  When  explictt  FD  method  is 
used,  which  is  the  most  laopular  FD  technique  to  date,  the  wavefield  of  a  specific  time  install  is  solved 
one  grid  point  by  one  grid  point  with  nearby  grid  data  at  previous  steps.  For  schemes  that  use  second- 
order  approximations  to  the  temporal  derivative,  only  two  grid  planes  of  rfisplacement  (or  stress,  velocity) 
must  be  stored  to  perform  the  updating  process.  Once  the  entire  grid  is  updated,  FD  then  proceeds  to 
compute  the  wavefield  of  the  next  time  instant  until  a  certain  preselected  number  of  iterat'nns  is 
reached.  The  output  of  FD  method  can  be  snaii^hots  of  the  entire  grid  at  spedTic  times  or  synthetic 
seismograms  recorded  at  specific  grid  points. 

‘'Fd2”  is  a  software  package  developed  at  Teledyne  Geotech  Alexandria  Laboratories  (TGAL)  dur¬ 
ing  the  past  several  years  for  generating  synthetic  seismograms  and  displaying  the  wavefiekts.  This 
package  consists  mainly  of  a  2-dimensional  2nd-order  explicit  linear  finite-difference  (LFD)  code.  In  this 
topical  report,  we  first  review  in  detail  the  algorithms  we  developed  for  the  topography  handling.  We  then 
present  sample  scripts  illustrating  the  usage  of  ‘id2"  as  well  as  several  supporting  routines  for  plotting 
out  the  synthetics,  generating  media,  and  the  graphic  visualization  of  wavefiekts.  Thus  this  topical 
report  serves  as  both  a  programmer’s  guide  and  the  user's  manual.  The  supporting  routines  included  in 
this  manual  are  ‘idjnenu".  'idspit",  "arrayplt",  “trpit",  “movie",  and  “fymedS",  which  constitute  the  most 
frequently  used  portion  of  our  software  package. 

Although  LFD  codes  are  becoming  very  popular,  our  code  does  posses  several  features  that  are 
not  that  standard.  Perhaps  the  most  distinguishable  feature  of  this  code  is  that  it  allows  simulations  with 
fairly  rough  topography.  Topography  can  focus  or  defocus  incident  body  waves  and  can  convert  body 
waves  to  surface  waves  and  vice  versa.  In  this  tutorial  we  therefore  include  an  overview  of  our  topogra¬ 
phy  algorithms  which  are  based  largely  on  our  published  work  (Jih  et  al. ,  1988).  Basically  the  algorithm 
we  developed  is  an  improved  version  of  llan  (1977).  On  the  inclined  free-surface,  the  vanishing  stress 
conditions  are  implemented  to  a  rotated  coordinate  system  parallel  to  the  inclined  boundary  as  previous 
works  did.  For  each  transition  point  on  the  topography  where  the  slope  changes,  we  use  the  first-order 
approximation  of  boundary  conditions  in  a  locally  rotated  coordinate  system  in  which  the  normal  axis 
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always  coincides  with  the  bisector  of  the  comer.  These  extrapolation  formulae  are  consistent  with  boun¬ 
dary  coixlitions  to  first-order  accuracy  in  spatial  increment,  the  same  as  the  classical  one-sided  explicit 
approximation  scheme  widely  used  for  the  flat  ftee-suiface  case.  Testing  resuHs  indicate  that  this 
scheme  works  stably  for  fairly  complicated  geometric  shapes  consisting  of  ridges  and  vaHeys  wtth  steep 
and  gentle  slopes  over  a  range  of  Poisson  ratios  of  practical  interest,  thus  enabing  us  to  study  more 
realistic  problems  for  which  the  topography  plays  a  significanl  role  in  shaping  the  wavefleld  and  analyti¬ 
cal  solution  might  not  be  available. 

A  few  words  on  the  evolving  history  of  Geotech's  LFD  code  are  necessary:  Z.  A.  Der,  J.  Bumetti, 
and  T.  McEHresh  initialized  the  prototype  code  design  in  late  seventies.  During  1978-1981,  they  imple¬ 
mented  the  2nd-order  explicit  FD  scheme  with  the  heterogeneous  formulation  of  KeNy  ef  al.  (1976),  the 
monochromatic  PISV  planar  source,  as  well  as  the  symmetric  boundary  condition.  K.  L.  McLaugNin,  T. 
McElfresh,  and  L.  Anderson  implemented  Ohnaka  (broadband)  PtSV  source,  the  point  (line)  source  and 
the  absorbkig  boundary  condition  (Clayton  and  Engquist,  1977;  Emennan  and  Stephen,  1983)  during 
1983-1985.  R.-S.  Jih  coded  up  the  fundamental-mode  Rayleigh  wave  generation  routine  (Boore,  1970; 
Munasinghe  and  Famell,  1973),  the  Ist-order  fonnuiation  of  free-surface  boundary  condition  to  handle 
the  polygonal  topography  (Jih  ef  al. ,  1988),  the  marching  grid  technique  for  extending  the  propagation 
distance  in  the  lateral  direction  (Jih  at  al.,  1989),  as  well  as  a  strain  filter  for  far-field  body  wave  synthet¬ 
ics  (Jih  at  al. ,  1989).  The  current  version  of  FD  code  ‘id2”  is  totally  different  from  all  earlier  versions  of 
Geotech's  LFD  codes:  ^'fdabct"  through  “fdabce")  after  a  series  of  heavy  revisions  by  R.-S.  Jih  during 
1986-1993,  even  though  several  subroutines  still  retain  their  original  names. 


2 


2.  TOPOGRAPHIC  ALGORITHMS 


2.1  Basic  Concepts  Underlying  Topographic  AlgorMims 

Consider  a  half-space  with  an  arbitrary  polygortai  free-surface.  DMne  a  Cartesian  coordinate  sys¬ 
tem  with  X  horizortai  (positive  to  the  right)  and  Z  vertical  (positive  downward)  as  in  Figure  1.  Assume 
that  the  material  is  perfectly  elastic,  isrrftopic,  and  homogeneous  with  oompressionai  and  shear  velodties 
a  and  p  respectively.  Let  U  and  W  be  the  X  and  Z  components  of  displacement.  Then  the  wave  propa¬ 
gation  at  interior  poirtfs  is  governed  by  the  following  hyperbolic  partied  differential  equations: 


^  ^2a^U(x.2.t)  ^  p2^  3^V(x.z.t)  ^  p2a^U(x.z.t) 

3t^  dx*  dxdz  3z^ 


or  ox^  oxoz 


..g3^(x.z.t) 


[lal 

[lb] 


A  grid  is  imposed  on  the  X-Z  plane  wtth  fixed  mesh  size  Ax  and  Az  .  We  let  x  ^  iAx,  z  =  kAz  ,  and 
t  =  lAt ,  where  At  is  the  increment  step  in  time,  and  i.  k.  and  I  are  positive  integers.  We  denote  by  Ujjg 
and  W|,ig  the  components  of  displacement  at  the  disaetized  grid  point  djo  ■  (lAx,  kAz,  lAt)  tfi  the  tth 
time  step.  By  replacing  the  partial  derivatives  in  equation  (1]  by  central  differences,  the  following  stan¬ 
dard  explicit  second-order  formulae  are  obtained: 

U|,k.kl  =  -  UlXH  +  2Ugg  +  (  Umjo  “  2U|,iu  +  Ui-ijg  ) 

+  ~  2Uixi  +  Uiji-i; ) 

+  -  P")  •  (  -  Wm -  Wm*.u  +  [1c] 

and 

Wijcj,,  =  -  W,^,  +  2Wu^  +  0t2(^)*(  W,j,^,j  -  2W,*j  +  W,j._,4 ) 

+  P"(^)*(  Wmxi  -  2W,J.J  +  Wmxi  ) 

+  •j~;^~^  (  “  P^  •  (  -  Ui+i,k-ij  -  +  UM.k-i.i)  •  [Id] 


The  reader  is  referred  to  ARerman  and  Kara!  (1968),  Boore  (1972),  and  Kelly  ef  a/.  (1976)  for  a 
more  detailed  discussion  of  the  approximation  of  partial  derivatives  by  finite  difference  discretizations. 

The  vanishing  of  stress  on  the  flat  free-surface  boundary  is  expressed  by  Xxz  =  0  and  Xzz  =  0  ,  or, 
equivalently, 
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{2a] 


^  3W 
dz  ^  dx 


=  0 


and 


-s£)iu 


=  0 


[2b] 


Alterman  and  Karal  (1968)  presented  an  explicit  central-dMerence  approximation  of  [2a,  b]  which  is 
stable  when  >  0.3: 

Uu  =  Uijt+2  -  -  Wm^i)  [2c] 

and 

Wu.  =  Wy^  -  ^  ( 1  -2-^)(U„,  .  (2dl 

Another  approach  which  adopts  one-sided  differences  was  proposed  by  Alterman  and  Rotenberg 

(1969); 

U.*  =  ^  (WwMi  -  W^>-,)  [26] 

and 

Wy.  *=  W,*,,  -  -^  ( 1  -2^)(U^,  -  U,.,j^,) .  t2fl 


At  the  vertical  surface  the  boundary  values  can  be  derived  by  the  following  transformations: 

u— >w.  w— ♦u.ax— »az.3z— »3x.  pj 


Although  the  one-sided  approximation  [2e.f]  is  appHcable  over  a  somewhat  more  restrictive  range 
of  Poisson’s  ratios.  some  cases  resuKs  are  even  better  than  the  centrakfifference  approximation 
[2c,d]  or  a  composed  second-order  scheme  by  Han  and  Loewenthal  (1976).  AH  algorithms  dtecussed  in 
later  sections  wHI  use  either  equations  [2c, d]  or  equations  [2e,f]. 

Alterman  and  Loewenthal  (1970)  proposed  two  approaches  to  the  calculation  of  the  displacements 
at  90^  and  270°  comers.  The  first  approach  recprires  that  the  normal  stress  components  on  both  the  hor¬ 
izontal  and  the  vertical  surface  be  zero,  and  in  place  of  condition  [2a],  it  imposes  a  more  restrictive  con¬ 
dition  of  dU/dx  -  dW/dz  «  0.  The  second  approach  srrxroths  the  comer  slightly,  making  the  tangent  to 
the  boundary  be  at  an  angle  of  45°  to  both  axes.  Two  similar  difference  schemes  following  the  second 
approach  can  be  found  Han  ef  a/.  (1975)  and  Fuyuki  and  Matsumoto  (1980).  Boore  et  al.  (1981) 
reported  that  Fuyuki  and  Matsumoto's  difference  equations  are  somewhat  better.  We  will  follow  the  idea 
of  the  second  approach  to  derive  formulations  for  arbitrary  comers. 

To  implement  arbitrary  topography  in  the  2-dimensionai  finite-d'ifference  scheme,  it  suffices  to 
approximate  the  topography  by  polygons.  Six  algorithms  are  presented  here  to  implement  the  free- 
surface  boundary  condition  for  the  six  separate  cases  shown  in  Figure  1  labeled  (A)  through  (F).  Only 
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the  formulae  for  increasing  elevation  with  increasirx)  x  are  given  for  brevity.  Since  the  calculation  of 
comer  points  always  requires  dispiacements  of  poirtfs  nearby,  it  is  necessary  to  follow  a  specific  order  in 
each  time  marching  step  to  solve  the  boundary  values:  first  comptrte  points  on  the  horizortal  and  verti¬ 
cal  free-surface  segments  with  (2c,..,  2f]  (or  transformation  3),  then  points  on  the  ramps,  and  finally  the 
comer  points.  The  topography  is  restricted  such  that  each  segment  of  the  polygon  consists  of  at  least 
three  points,  so  that  comer  points  are  always  separated  by  ramp  points  (at  least  one).  The  selection  of 
spatial  sampling  increment.  Ax  and  Az,  depends  on  the  velocities,  frequencies,  and  geometric  resolution 
of  the  problem  to  be  solved.  For  complicated  topography  it  is  necessary  to  adopt  smaller  Az  (relative  to 
Ax)  to  ensure  the  inclusion  of  the  most  gentle  rarr^  encountered,  namely  the  ramp  with  slope  Az/Ax. 
The  temporal  spacing  At  is  then  determined  by  the  mesh  spadngs  as  well  as  by  wave  velocities  so  as 
to  match  the  general  stability  constraint. 

The  underlying  procedure  for  all  the  six  topography  algorithms  is  the  following: 

(1)  At  every  node  on  the  topography,  set  up  a  rotated  coordinate  system  X-Z'  which  accounts  for  the 
local  geometric  orientation.  For  comer  points,  simply  take  the  bisector  of  the  comer  as  the  Z'  axis. 
For  points  on  the  ramp,  this  means  that  Z'  is  the  local  normal  to  the  ramp. 

(2)  Set  up  a  four-atom  molecule  with  diagonals  parallel  to  the  X-Z'  axes  and  let  the  unknown  node 
<i,k>  be  one  of  the  atoms  by  default.  Remaining  atoms  are  selected  so  that  their  displacement 
can  be  computed  with  interior  points,  thus  ensuring  the  scheme  remains  explicit.  The  molecule 
itself  should  be  as  close  as  possible  to  the  boundary  to  reduce  error  due  to  coarse  sampling  of  the 
wave  field.  Usually  it  is  convenient  to  put  an  atom  just  one  grid  below  node  <i,k>. 

(3)  Apply  equations  [2c,d]  or  equations  [2e,f]  to  the  four-atom  computational  molecule  to  solve  for  the 
boundary  value  at  node  <i,k>  in  the  X'-Z'  cooRiinate  system.  The  solution  is  then  rotated  back  to 
the  original  X-Z  system. 

The  details  of  the  various  algorithms  are  presented  in  Section  2.3. 
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2.2  Han’S  Method  for  Constant-slope  Ramp 

Man  (1977)  uses  a  rectangular  grid  with  fixed  height  Az  and  variable  width  Axj .  Each  AXj  is  deter¬ 
mined  via  Ax  ■  Azcote  ,  where  6  is  the  iocai  angle  between  the  X-direction  and  the  inclined  tree-surface. 

Suppose  that  node  <i,k>  lies  on  an  irx:iined  tree-surface  with  slope  Az/Ax  (Figure  2).  It 
0  ■  tan~VA2/Ax)S45"  .  Let 

U',>=U'o-sin20(W'M.k-W',j„.,)  I4a] 

and 

ViTij.  =W'o  -  Sin20(1  -  2^)(UV,.k  -  U^,)  .  [4b] 

where 

Uo  =  CqUmmi  +  ( 1  -  Eo  )U|j«+2  .  I5a] 

Wo  =  eoWMx.i  +  (1-eo)W,j^.  [Sb) 

and 

Co  ■  2sln^  .  [5cl 


If  0  >  45** ,  then  replace  equations  [5]  by 


Uo  =  CqUmmi  +  ( "•  “  Eo  )Umjj  . 

(6al 

Wo  =  eoWi^tjH-l  +  (  1  “  Eo  > 

[6bl 

Co  =  2cos^ . 

[6cl 

Finally  rotate  U'j,K,  W'i,k  of  equations  [4]  back  to  the  original  system. 

We  have  observed  that  without  appropriate  manipulation  of  the  transition  points  on  the  topography 
as  proposed  in  the  present  work,  this  algorithm  alone  would  suffer  from  instabilities.  Furthermore,  this 
single  algorithm  with  varying  grid  spacing  approach  yields  more  complexity  in  its  use  than  does  the 
multi-algorithm  with  fixed  grid  spacing  described  in  the  present  work. 
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2.3  Algoilthms  Us«d  in  FD2 
(A)  Constant  Steep  Slope 


Suppose  that  a  grid  point  <i.lc>  lee  on  an  incfined  free-suiface  with  slope  m&z/Ax,  mi 2  (Figure 
3).  H  we  rotate  the  X-Z  coordinate  ^em  by  the  anj^e  e  ■  tan~'(m  az/ax)  oounterdockwise  to  the 
X -Z'  system,  then  in  Hght  of  equations  (2a.b]  we  have 


dir 

dz' 


+ 


[A-la] 


and 


3W' 

3z' 


+  (1-2 


^  dUl 
o^'ax' 


=  0. 


[A-lb] 


where  U'.  W'  are  the  dispiaoement  in  X'.  Z'  direction  respectively-  These  two  equations  wW  be  used 
to  solve  for  the  displacement  at  every  grid  poM  on  the  free-surface  topography  with  different  0. 


Let  point  0  be  the  projection  of  <i.k>  on  the  line  Joining  points  <i+1  ,k>  and  <i  A+m>.  With  first-order 
acojracy,  the  displacements  at  0  can  be  approximated  by  the  folowing  linear  interpolations: 

Uo  =  eoUM,k  +  { 1  -  So  )Uuh<ii 

and 


Wo  ^  “  E©  )W|j|*«n  I 


where 


Co  ■  Sin^ . 


Rotate  the  displacements  at  the  three  atoms  0,  <i+1.k-1>.  and  <i,k+m'1>  to  the  inclined  system 
X-Z'  and  substitute  them  into  equations  [A-l].  We  have  the  following  extrapolation  formulae: 

U',*  =  ITo  -  cos0sine(W'M*-,-W',j^,)  IA-2al 

and 

W'.j,  =  W'o  -  ooseslne(1-2^)(U'Mj^i-U',^,) .  [A-2b] 


U'i>  and  W'i  k  are  then  rotated  back  to  the  original  X-Z  system  by  the  angle  -0  . 
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Figure  3 


(B)  Constam  GMMIa  Slope 

For  grid  point  <iJo  on  an  inclined  free-surface  with  slope  Az/Ax  (Figure  4),  let 

8  ■  tan~'(AzrAx) 
and 

eo>8in^ . 


Let  point  0  be  the  projection  of  <i.io  on  the  line  joining  poMs  <ii-1  M>  and  <ij(-fl>.  Again,  the  displace- 
mer4s  at  0  can  be  approximated  via 

Uo  =  eoUi^t  +  ( 1  -  Co  )Uiji4.i 


and 


Wo  =  1  )W|j^i  . 


Now  in  the  X'-Z'  coordinate  system,  displaoements  t/iji ,  W'lj,  can  be  solved  via  equations  [B-1]  as 

U'u  =  ITo  -  sinecose(W'w.k  -  W'u,,)  [B-la] 

and 

W'lj.  -  W'o  -  sineoose(1  -  2^)(UV,j,  -  .  [B-lb] 

We  then  rotate  l/u,.  W'l^c  back  to  the  original  X-Z  system  with  the  angle  -0.  Note  that  this  algo¬ 
rithm  is  simpler  than  llan’s  (1977)  formulation,  and  the  weighting  factor  sinecosO  in  formulae  [B-1]  is 
applicable  whether  6  ^  4S**  or  not.  although  in  most  actual  applicalions.  8  is  rarely  larger  than  45”.  Also 
note  that  we  use  gentle  and  steep  in  algorithms  A  and  B  simply  to  distinguish  whether  m,  the  eievation 
rising  from  grid  column  i  to  column  h-l  in  terms  of  Az,  is  equal  to  one  or  greater  than  one. 
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(C)  Concav*  HorizonttiM<>a«ntl«  Slope  Transition 

For  grid  points  at  the  bottom  of  a  canyon  wtoi  slope  mOz/Ax  .  m  2 1  (Figure  5).  assuming  tttat 
Ax  2;  Az  ,  let  e  ■  tan~’(mAz/Ax),  eo  ■  tan(a^)Az/Ax,  e,  ■  tan(e/2)co«e.  and 

Uo  =  (1-eo)Uu»i  +  jw-1  •  ®1 

Wo  =  (1-eo)W|j,*,  +  SoWkijm-i  .  (C-lbJ 

U,  =  (1-e,)UUij.  +  e,Un  .  (C-lc] 

W,  =  (1-€,)Wm*  +  .  lC-1dl 

U2  =  ('l-*i)Uk-i*  +  SlUwijMii .  IC-le] 

and 

Wa  =  (l-ei)Wwj.  +  •  IC-lfl 

Then, 

U'ij.  =  U'o-^(W'2-W',)  lC-2a] 

and 

w'l*  »  w'o  -  ^(1  -  -  ir,) .  (C-2bi 

Note  that  we  rotate  U,  W  by  the  angle  W.  instead  of  0  to  get  IT.  W'  so  that  the  T  uds  is  con¬ 
sistent  with  the  directicn  of  the  local  bisector  at  grid  point  <i,k>. 

For  the  case  Ax  <  Az  ,  substitute  co  ■  cot(6r2)Ax/Az  and  Ci  ■  tane  tan(e/2)  in  formulae  [C-1], 
and  replace  Az/2Ax  to  equations  [C-2]  by  oot(0/2)/2  . 
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Figure  5 


(D)  Concav*  G«fitl«-to>«teep  Slop*  Transition 


For  grid  points  with  iefi  slope  nAz/Ax  and  rigM  slope  mAz/Ax,  m  >  n  ,  assuming  that  Ax  ^  nAz  ,  let 
e  m  tan~'(nAz/Ax),  ^  ■  tan~'(mAz/Ax),  i)  « (6  +  ^  )/2.  cq  ■  ei  ■  tanOcohi.  62  >  tarn)  cot#,  and 

Uo  =  (1-Co)U»n*  +  eoUw-ijMfl  . 

(D-1  a] 

Wo  =  {1-€o)Wv4.^ji  +  1 

(D-1bl 

Ui  =  (l-Ci)UuMn  +  CiUl-ijHn  . 

(D-lcl 

W,  =  (1-e,)WMj««  +  e,WMj««. 

(D-1dl 

U2  ~  (1— CiUkl.lnn  > 

(D-1  el 

and 

W2  =  (1-ei)W4.iji  +  StW^XHn  • 

(D-in 

Then, 

u'l  >  =  u'o  -  7 — V-r(W'2  -  w'l) 
tanri+tane 

(D-2al 

and 

W',  w  =  W'o  -  7 — ~r(1  -  2-^)(U'2  -  U'i) . 

'•*  "  tami+tane  0* 

(D-2b) 

Note  that  we  rotate  U,  W  by  angle  i\  to  get  U',  Vyr  so  that  the  Z'  axis  is  consistent  with  the  direc¬ 
tion  of  the  local  bisector  at  the  grid  point  <i,k>. 

For  the  case  that  Ax  >  nAz  (Figure  6),  substitute  Cq  ■  tanetamt,  ■  tanOcotn.  and  tz  m 
in  formulae  [D-1],  and  replace  tanTi+  tane  in  [0-2]  by  coht4COte.  Also,  replace  [D-ia,b]  with 

tamtcot#. 

Uo  —  (^~^)Ui4M«  +  EoUmJho  > 

(D-3a] 

Wq  =  (1~eo)Wj_|(^  +  eoW4.iji4fl  . 

[D-3bl 
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(E)  Convex  Change  In  Sope 

For  grid  points  with  decrease  in  slope  from  nAz/Ax  to  mAz/Ax,  n > m,  let  Ob tan~’(nAzfAx), 
4  m  tan'^imAz/Ax),  y\m{Q  +  ^y2,  and  sq  ■  taraiAz/Ax.  H  cq  is  less  than  1,  then  let 

Uo  =  (1-«o)U|j,+t  +  CqUu.ijm'I  [E-fa] 

and 

Wo  *  +  eoWu.iX4-i  . 

as  shown  in  Figure  7b  where  atom  0  Hes  between  <i.k-i-1>  and  <i-i-1,k-f1>. 

Uo  =  (1-e3)Ufijin  + 
and 

Wo  =  (1-e3)WMMi  +  esWMj. .  lE-ld) 

so  that  atom  0  lies  between  <1+1  ,k+l>  and  <i+l.k>  (Figure  7a). 

Now  depending  on  the  value  of  m .  we  have  2  cases  : 

(E.1)  m  =  1.  Let  <i+1.lo  be  an  atom  (Figure  7a).  Then  the  riagonal  of  the  molecule  passing 
through  <i+1,k>  wiH  intersect  the  ith  grid  column  somewhere  between  <i,k>  and  <i.k+n>.  Let  the  inter¬ 
section  point  be  atom  1  (Figure  7a).  Note  that  tan#  <  taraf  <  tanO,  hence.  ei  ■  (tamicot#  -l)/(n-1)  is 
always  between  0  and  1 .  The  displacement  at  atom  1  can  be  interpolated  as 

Ui  =  (1-et)Ugw.i  +  SiUijMii  IE-2al 

and 

Wi  =  (1-c,)W,j,*,  +  eiW,j^^.  (E-2bl 

Note  that  we  use  point  <i,k+1>  rather  than  <i,k>  in  equations  (E-2),  since  the  displacement  at  <i,k>  is 
unknown.  Then, 

U'u«=U'o-e2(W'MjrW'i)  IE-3aJ 

and 

W'ij,  =  W'o  -  ezd  -  2^)  .(UVix-U'i) .  tE-3bl 


{E-1bl 

Othenwtee  let  ea  b  i/co,  and 
[E-lcl 


where  ez  s  tan#  if  eo  ^  1.  or  coin  if  Co  >  1- 

(E.2)  m  >  1.  Let  atom  2  be  at  node  <i,k+1>  (Figure  7b),  and  a  (2Azcotn/Ax),  so  ei  is  always 
between  0  and  1. 


Ui  =  (1-ei)Ujj,_i  +  e,Uj+ij,_i 


IE-2C] 


and 


W}  =  (1-Cj)W’j|,_i  +  e|Wf,.ij,_j . 


IE-2d] 
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The  displacetneiUs  at  poim  <j,k-1>  are  computed  by  extrapolation  (see  algorithm  (G))  since  rs  out¬ 
side  the  grid.  Now, 

ir.*  =  l/o  -  eaCVir,  -  IE-3c] 

arKl 

Vir,*  =  ViTo  -  e2(1  -  2^){U\  -  U'i*.,)  .  IE-3dl 

a 

where  ca  ■  O.Stami  if  Eq  ^  1.  or  ea  ■  Ax/(2Az)  if  eo  >  1  . 

Note  that  again  we  rotate  U,  W  by  angle  i)  to  get  (/.  VIT  ,  so  that  the  Z'  axis  is  consistent  with  the 
direction  of  the  local  bisector  at  the  grid  point  <i.k>. 
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(F)  Convex  Gentle  Siopo4o*horlzontal  Transition 


Without  loss  of  generality,  we  may  assume  that  the  topography  changes  slope  from  Az/Ax  to  0  at 
the  crests  (Figure  8).  Let  ^  ■  tan~'(Az/Ax),  cq  ■  tan{^2)tan^,  and  ei  >  tan(4/2)cot^.  In  this  case  the 
offset  between  the  Z-direction  and  the  Z'-direction  is  ^2.  Let  the  intersection  of  the  local  normal  wtth  the 
(k-i-l)-th  grid  row  be  atom  0.  The  other  two  atoms  are  chosen  so  that  they  are  collinear  with  node 
<i,k+1>.  Note  that  eq  ^  1  provided  Az  ^  ^  Ax  .  /.e. ,  ^  ^  60°.  Let 


Uo  =  (1-€o)Lli4t+i  +  eoUH-i,iM-i  . 

IF-1a] 

Wo  =  (1-€o)W|j,+i  +  EoWh.iji+1  , 

[F-lb] 

Ui  =  (1-Ei)Ui_i>+1  +  EiUi_i.k+2  1 

(F-lcl 

W,  =  (1-ei)WMMl  +  £lWM>.2. 

lF-1d] 

U2  =  (1~ei)LlH.ij(+i  +  £iLli+i,k  1 

lF-1el 

W2  =  (1-€t)Wi4.iji+i  +  EiWi+i.k  • 

[F-in 

Then, 

U'u  =  U'o  -  ^(W'2  -  W',)  [F-2a] 

and 

=  W'o  -  -^{1  -  2^)(U'2  -  U',)  .  IF-2b] 

Here  we  rotate  U,  W  by  angle  ^2  to  get  U',  W'  so  that  the  Z'  axis  is  again  consistent  with 
the  direction  of  the  local  bisector  at  the  grid  point  <i,k>. 
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(G)  Fictitious  Boundary  Points  for  Steep  Slope 


Suppose  that  a  grid  point  <j,k>  lies  on  the  tree-surface  with  slope  mAz/Ax,  mt  1  .  Due  to  the 
explicit  scheme  (Equations  [1c]  and  [id])  we  adopt  for  the  interior  points,  the  calculation  of  the  dKspiaoe- 

ments  at  inner  points  <1+1. k>.  <l+1,k-1> . <i+1J(-m+1>  requires  displacement  values  at  <IJt-1>.  <IJt- 

2>,...<i,k-m>  which  are  outside  the  topography.  This  dfficuky  can  be  easily  removed  by  evakiating  m 
fictitious  layers  above  the  topography  as  follows  (Figure  9): 

For  1  ^  n  s:  m,  let  eo  ■  (n/m)sin2e.  where  6  ■  tan-’mAz/Ax,  use  Uo  =  eoUMji-m  +  (1  -  eo)Uy,  and 
Wo  -  eoWMj(-m  -t-  (1  -  eo)W|^  as  the  approximate  displacement  of  the  orthogonal  projection  of  point  <i,k- 
m+n>  on  the  free-surface. 

Now  rotate  displacements  at  points  <i+l  ,k-m>,  <i,k>  and  0  by  angle  0  as  before  to  extrapolate  the 
approximate  motion  at  point  <i,k-m+n>  : 

=  u'o  -  •^sinecose(W'Mji-n.  -  w'ljj 
and 

=  W'o  -  (l-2^)^slnecose(UVaHn  -  . 

Then  rotate  U',  W'  back  to  the  original  system. 
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2.4  Some  Ranuufcs  on  Topographic  Algorithms 

Even  in  the  simplest  haV-space  models,  the  central-diference  approximation  (foimulae  2c, d)  and 
the  one-sided  approximation  (formulae  2e,f)  of  the  boundary  conditions  are  only  conditionally  stable. 
When  quarter  planes  are  involved,  the  range  of  stability  is  inevitably  smaller.  Since  our  schemes  are  a 
(rotated)  extension  of  these  two  methods  and  since  we  are  deaBng  with  more  compHcated  irregular 
geometries,  a  stronger  limitation  of  Poisson's  ratios  should  apply,  it  is  obvious  that  our  scheme  has  trun¬ 
cation  error  of  order  one  (with  respect  to  the  mesh  spadngs)  due  to  the  spatial  disoetizalion  and  linear 
interpolation,  as  well  as  the  local  asymmetry  in  our  formulation,  even  though  the  original  (f.e. ,  unrotated) 
central-difference  approximation  (2c,d)  has  secorKf-oider  accuracy.  Nevertheless,  better  accuracy  can 
be  achieved  by  sampling  the  polygon  with  a  finer  grid  at  the  expense  of  computation.  All  algorithms  are 
of  the  explicit  type,  since  only  interior  nodes  nearly  are  needed  for  evaluaiirtg  tangential  derivatives. 
The  actual  surface  configuration  implied  by  this  scheme  is  a  somewhat  smoothed  version  the  original 
discretized  polygonal  boundary. 

Recently  a  number  of  techniques  have  been  pursued  in  an  effort  to  improve  the  computational  per¬ 
formance  of  finite-difference  solutions  to  wave  equations.  These  include  higher-order  schemes,  implicil 
rather  than  explicit  methods,  velocity-stress  schemes,  etc..  Our  fonnulations  were  developed  mainly  to 
adapt  to  the  traditional  explicit  second-order  finite-dtiference  schemes  for  wave  equation  (equations  1) 
as  described  by  Alterman,  Boore,  and  Kelly  because  of  its  popularity.  Incorporation  of  similar  treatments 
of  a  free  boundary  with  elasfodynamic  equations  might  be  promising,  since  a  velocity-stress  scheme 
{e.g. ,  Virieux,  1986)  is  more  tractable  than  the  standard  cfisplacement  schemes  at  liquid-solid  interfaces. 

ANemative  approaches  are  available  to  handle  seismologicai  problems  involving  topographic  stroc- 
tures.  For  instance,  the  finite-element  method  offers  another  natural  solution  with  which  one  can  use 
more  general  grid  elements  such  as  triangles  or  even  higher  order  forms  instead  of  the  rectangular  grids 
used  in  the  finite-difference  method.  Drake  (1972)  and  Ohtsuki  et  al.  (1984)  use  the  steady-state 
finite-element  method  and  the  transient  finite-element  method  respectively  to  study  the  scattering  of 
Rayleigh  waves  by  topography.  Finite-difference  and  finite-element  methods  can  be  interpreted  as 
degenerate  cases  of  each  other  under  some  conditions  (Lapidus  and  Finder,  1982).  It  is  reported  that 
finite-element  methods  give  superior  accuracy  relative  to  finite-difference  methods  when  modeling  elastic 
media  with  Poisson's  ratio  between  0.3  and  0.45,  and  the  accuracy  is  comparable  to  finite-difference 
methods  for  media  with  Poisson's  ratio  less  than  0.3  (Marturt,  1984).  The  major  advantage  of  the  finite- 
difference  method  over  finite-element  methods  is  its  simplicity  in  implementation.  Our  work  shows  how 
the  finite-difference  technique  can  be  extended  to  an  irregula^  free-surface  with  special  treatments. 


3.  SAMPLE  INPUT  FILES  FOR  FD2 

“Fd2'‘  reads  in  control  parameters  ttwough  the  standard  input  and  dwnps  the  error  messaoM  to 
the  standard  output.  The  output  synthetics  as  %wel  as  the  snapshots  of  the  wavefieU  are  stored  in 
separate  output  buffers.  The  source  medLim  file  is  specified  in  the  input  file  or  assumed  to  be  a  preset 
homogeneous  haK  space.  Topography  can  be  specified  as  an  input  file.  A  sample  input  file  is  shown 
here  with  each  line  of  argument  cfiscussed  in  more  detail  at  the  end.  FCr  the  user's  convenience,  an 
interactive  routine  "Ujnenu"  is  avalable  which  woifid  create  the  input  parameter  file. 

Sample  Input  Parameter  Hie 

[I]  grid  dimenaion  km,kh 
800600 

P]  K2  tfiadnga  at  ihm  grid  math  A  tamfioni  apadng:  dx.di.dl 
0.25  0.250.(05 

PJ  homoganaoua  P).  hataroganaeus  (1.2i 
2 

inVp 

biVa 

InRho 

fdjtopognphym 

none 

15]  ehoote  wave  type 
3 

(6]  incident  mgle,  iOjO,  wvins 
20.00400  325.0 

[7]  option  for  enapehot  dieplay:  oompanentAOC.begin  8  end  column 

re-  1 

1  800  1800132200 

401  800  1800132  100 

201  00  15001  2600 

301  00  801  2000  2250 
701  00  801  1261  2  100 
901  00  1201  1561  2  100 
1101  00  1381  1801  2  100 
1301  00  1681  2161  2  100 
1501  00  1941  2341  2  100 
P]a  of  aeiamographa 
10 

P]  coordmatea  of  senaors:  (2  tinea,  tne  I^X's,  tine  2mZ's) 

1  200  400  600  800 -1 ‘200 -400 -600 -800 

2  2  2  2  2  3  3  3  3  3 

[10]  total  a  of  derations:  *  of  deradons/snapshot 
8001  200 

[II]  marching  grid  (>Om^yes)  ?  total  grid  widd>  deskad: 

400  1200 

Remarks: 

Item  (1) 

The  grid  size  specifies  the  number  of  nodes  in  the  horizontal  and  vertical  directions.  The  program 
will  stop  if  the  requested  grid  size  exceeds  the  predefined  size,  which  is  currently  compiled  with 
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800  by  600  nodes  on  CSS*  SUN4. 

Item  (2) 

The  grid  mesh,  AX  and  AZ.  is  measured  in  km.  The  temporal  spacing  is  in  sec.  Once  the  grid 
mesh  is  specified,  the  temporal  qaadng  should  be  determined  so  that  the  stablily  can  be  ensured 
throughout  the  iterations, 
item  (3) 

if  0.  the  default  homogeneous  model  with  density,  p.  of  2.5gnVcc.  and  X  «  p  -  21  for  the  soil, 
which  implies  Vp  -  5.02km/a  and  Vs  •  2.9knVs. 

If  1.  the  program  expects  to  read  in  X.  p.  and  p  fields. 

If  2.  the  program  expects  to  read  in  Vp.  Vs.  mxf  p  fields. 

Item  (5) 

The  initiai  wave  could  be 

broadband  planar  P  wave  of  Ohnaka  shape, 

(•1)  broadband  planar  S  wave  of  Ohnaka  shape. 

(•4-2)  monochromatic  planar  P  wave  of  sinusoidai  shape. 

(•2)  monochromatic  planar  S  wave  of  sinusoidal  shape, 

(43)  pure  oompressional  wave  generated  at  a  single  point, 

(-3)  double-oouple  point  source. 

(44)  fundamertal  mode  Rayleigh  wave  w^  Ricker  wavelet  shape, 

(45)  a  (time'Series)  P  driver  file  shaking  a  single  point, 

(46)  arbitrary  wavefieid  setting, 

(4?)  broadbarKi  planar  P  wave  of  Gaussian  shape, 

(•7)  broadband  planar  S  wave  of  Gaussian  shape, 

(48)  ripple-fired  explosions  with  multiple  sources  of  type  (43), 

(49)  ripple-fired  explosions  wtth  multiple  sources  of  type  (45). 

Except  for  cases  (45)  arxf  (49)  in  wNch  the  shaker  or  driver  ffle  must  be  generated  elsewhere  in 
advance,  “fd2”  is  completely  self-contained  to  initialize  the  wavefields  at  2  consecutive  time 
instarfis  to  produce  proper  wave  propagation. 

Item  (6) 

(a)  the  wave  length  (in  km)  is  needed  for  a  point  source,  sine  or  Rayleigh  wave 

(b)  the  location  of  the  source  point  can  also  be  used  to  ac^ist  boundary  conditions 

Absorbing  boundary  condition  is  default  for  the  bottom  and  side  edges.  For  rwrmally  incidert 
planar  waves  (options  1,  2,  7),  the  side  edges  would  be  symmetric.  Free  surface  is  assumed  on 
the  top  of  the  grid  whenever  a  nontrivial  topography  is  involved.  AM  these  boundary  conditions  can 
be  altered  by  choosing  appropriate  input  parameters.  For  instance,  in  the  case  of  pokit  source 
(options  3,  5,  8,  9)  without  topography,  there  are  3  more  options  avaiiable  by  playing  with  incident 
angle: 

(1)  0-degree  indicates  that  all  4  edges  are  symmetric, 

(2)  360-degree  indicates  that  an  4  edges  are  absorbing, 

(3)  720-degree  indicates  that  a  symmetric  top  with  absorbing  skies  and  bottom. 

These  options  are  meant  only  to  demonstrate  the  effects  of  miscellaneous  boundary  conditions.  In 
reality,  models  with  free-surface  boundary  condition  on  the  top  and  the  absorbing  boundary 
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condition  at  the  remainino  three  edoes  would  be  the  most  usefui  ones,  in  thts  sample  input  Me, 
the  incidert  angle  of  20**  has  no  effect  and  hence  the  defauR  boundary  ooncMion  for  source  option 
•«-3  win  be  used,  which  is  absorbing  side. 

Item  (7) 

Program  *111(12"  converts  the  numetfcai  wavefields  into  the  character  wavefieidS  and  stores  these 
snapshots  in  ASCii  text  Mes.  The  output  wav^ieid  ooukJ  be  the  whole  grid  or  only  a  portion  of  the 
grid  as  specified  in  the  input  parameter  file.  An  toput  parameter  determines  whether  fixed  gain  or 
automaticaly  ati^ustaUe  scale  is  used  in  the  conversion  procedure.  Displaoements  andtor  strain 
may  be  recorded  as  a  time  aeries  at  any  interior  points  for  the  strain  or  any  grid  points  for  the  dis¬ 
placements. 

7  options  for  snapshot  displaying  are  available: 

H  :  horizontal  wave  fields 
V  :  vertical  wave  fields 
B  :  both  horizontal  &  verticai  wve  fields 
S  :  diiatational  &  rotatiortal  strain  fields 
M  :  combined  (fispiacement  fields 
A  :  an  components  wRI  be  shown 
AGC  >  0  ->  resolution  adjusted  for  each  frame  individually. 

At  most  9  slices  (stripes)  can  be  saved  in  compressed  files  with  prefix  “SI"  through  “S9",  respec¬ 
tively.  Each  stripe  needs  one  fine  to  specify  (1)  beginning  column,  (2)  endtog  column,  (3)  begin¬ 
ning  time,  (4)  ending  time,  (5)  beginning  row,  and  (6)  endtog  row.  If  any  of  the  three  endtog  entity 
is  smaller  than  its  associated  beginning  entRy,  then  that  stripe  wH  not  be  shown  at  aR.  in  this 
example,  only  two  stripes  are  stored  in  compress  ASCII  files.  Fles  SI*  store  the  upper  200  rows 
of  the  whole  wavefiekf  (800-grid  wide)  at  a  rate  of  one  snapshot  per  200  Rerations  {cf.  Rem  [10D 
throughout  toe  8001  Rerations  (cf.  Rem  (lOj).  FRes  S2*.  however,  store  the  cerRrai  portion  of  the 
grid  from  column  401  to  column  800.  Also,  only  the  top  100  rows  wW  be  shown,  as  specRied  in  the 
input  file.  The  remaining  seven  stripes  aN  have  an  ending  column  (00  in  this  case)  smaHer  than 
the  beginning  column,  which  simpiy  turns  (rtf  their  srrapshot  output.  AH  snapshots  are  shown  from 
row  2  and  down.  However,  a  dRfererR  “begirming"  row  can  be  given  to  force  the  gray  level  be  nor¬ 
malized  by  the  peak  displaoement  of  a  subgrid.  In  this  example,  the  gray  scale  is  determined  by 
32th  row  and  below,  aRhough  the  top  portion  wilt  be  shown  as  wen. 

Item  (9) 

(a)  list  X-coordinates  of  all  sensors  first 

(b)  negative  X-ooordinate  rrreans  strain  sensor 

(c)  R  should  be  consistent  wRh  toe  number  of  sensors  as  specified  in  Rem  (8) 

Item  (10) 

There  will  be  8001  Rerations  in  total  in  this  exarrple,  and  one  snapshrx  win  be  produced  for  every 
200  Rerations. 

Item  (11) 

If  a  nunrber  larger  than  0  is  given,  the  marching-grid  feature  will  be  turned  on,  and  the  grid  will  be 
shifted  by  that  number  of  columns  when  the  wavefront  is  approaching  the  right  edge.  In  this 
example,  400  columns  will  be  shifted  each  time.  The  total  grid  width  desired  is  1200,  which  is  the 
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original  grid  size,  800,  plus  400  columns,  so  only  one  shRing  wM  be  pertonnsd  (cf.  Step  4  of 
Example  3). 


Sample  Topography  File 

The  following  sample  file  defines  a  simple  ramp: 

10 

00  00  00  00  00  00  Of  02  00  <M  05  os  07  os  os  os  flS  os  os  os  M  os  os  M 

The  1st  Hne  gives  the  “referenoe  floor  level  of  the  topography.  The  2nd  me  gives  the  elevation 
with  respect  to  the  reference  level.  Fbr  example,  the  *1)3'*  in  line  2  actualy  means  the  7th  (alO-S)  row 
from  the  top  of  the  grid.  Note  that  the  reference  level  Ksel  must  be  deep  enough  so  that  the  peak  of 
the  topography  is  stM  within  the  grid.  Each  segment  must  contain  at  least  2  sub-segments  before  the 
slope  changes. 


Sampla  input  File  of  Valoeliy  llodel 


gridglim:  20  2S 
aMafiadno:  0.1000  0.1000 
aaKSMkr  kladkm 
Exauctu!  tton  uno(h$f  /mdW 
S  (dummy  tim) 

$  (dummy  Imf 
/(dummy  tm} 

$  (dummy  tnu) 

040276*01  0.40026*01  OM1006*01  081226*01  0.81206*01 
0807K*01  082136*01  0.82006*01  082206*01  0.82006*01 
088046*01  086086*01  088036*01  088426*01  082806*01 


Lines  1  and  2  specify  the  grid  size  and  and  the  ^  spacing.  Lines  3  and  4  are  lor  kfendflcatton 
purpose.  Lines  5  through  8  are  dummies.  The  remaining  Rnes  give  P-waye  velocity  (or  S-wave  velo¬ 
city,  or  density)  at  the  node  (i.j)  where  i«l,...,kw;  |-l,...J(h.  Lines  1  and  2  are  read  with  fixed  format  of  ( 
i1  Ox,  2i5 )  and  ( I3x,  2f  10.4 ).  respectively. 
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4.  OUTPUT  PILES  GENERATED  BY  F02 

A  number  of  different  output  formats  are  corrmonly  used  in  finle-ditference  calculations.  In  fact, 
one  of  the  advantages  of  the  finite-dWerence  method  is  the  ability  to  save  paflide  dtepiacements  for  any 
number  of  points  on  the  grid  and  for  al  times.  The  symhetics  are  nothing  but  the  time  history  recorded 
at  specific  grid  points,  for  example.  A  particularly  informative  way  to  visualize  wave  propagation  through 
the  models  is  by  using  numerical  Schlieren  diagrams.  These  “snapshots”  of  the  wavefieid  traveling 
through  the  grid  are  generated  by  saving  the  di^jiacement  at  each  grid  point  for  a  given  time  step. 

In  our  implementation,  the  program  converts  numerical  wavefields  irSo  character  wavefields  and 
stores  these  snapshots  in  ASCII  text  files,  it  consumes  relatively  smal  storage  space,  as  each  grid 
point  only  needs  one  byte.  The  output  wavefieid  could  be  the  whole  grid  or  only  certain  portion  of  the 
grid.  An  input  parameter  determines  whether  a  fixed  gain  or  an  automaticaly  adjustable  scale  is  used  to 
the  conversion  procedure. 

Displacements  and/or  strain  may  be  recorded  as  a  time  series  at  any  interior  points  for  the  strain 
or  any  grid  points  for  the  displacements. 

Program  also  stores  the  wavefiekts  at  two  consecutive  instarls  and  ail  required  parameters  at  a 
prespecified  rate  so  that  it  can  be  re-started  to  case  the  job  is  terminated  to  the  middle. 


5.  EXAMPLE  1:  PLANAR  WAVE  INaOENT  UPON  SIMPLE  RAMP 


The  first  example  demonstrates  how  to  set  an  experiment  for  modeSng  the  propa^tion  of  a 
normally  inddent  plane  P  wave  through  a  rmdel  with  a  45**  ramp  on  the  top  d  the  grid  and  von  Neu¬ 
mann  (/.e.  0-8tope  or  synwnetric)  boundary  condtions  on  both  sides.  The  homogerreous  medium  has 
compressionai  and  shear  velocities  of  5.02  and  2.898  km/s,  respectively. 


Step  1:  Generating  Topography  Hie 

car«  />  ktTOP 
66 

00000000000000000000000000000000000000000000 
00000000000000000000000000000000000000000000 
00000000000000000000000000000000000000000000 
00000000000000000000000000000000000000000000 
00000000000000000000000000000000000000000000 
00000000000000000000000000000000000000000000 
00000000000000000000000000000000000000000000 
00000000000000000000  01  0203040606  070600  10  11  12 
21  22  2324  2526272620  30  31  32  33  34  3636  S7  36  30  40  41  42 
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525354 

555657565060 

61 

6263 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 
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64 
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64 
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64 

64 

64 

64 

64 
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64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

e. 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 
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64 

64 

64 

64 

64 

64 

64 

64 
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64 
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00 

00 

00 

00 

00 

00 

00 

00 

00 

000000 

000000 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

0000 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

00 

0000 

00000000 

0000 

13 

14 

15 

16 

17 

18 

10 

20 

43 

u 

45 

46 

4748  40 

50 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 
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64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

64 

Step  2:  Running  FD 

#»-  generate  input  parameter  lie  for  U 
cat«  EOF>  inFD 

1)  grid  dimension  kw.kh 
500202 

2)  spadngs  of  the  grid  mesh  6  temporal  spacing:  dx,dz,dt 
0.01  0.01  0.001 

3)  homogeneous  (0),  hetarogeneaus  (1,2) 

0 

inVp 
in  Vs 
inRho 

4)  topography  fie 
inTOP 

5)  choose  wave  type 
7 

6)  alpha,  iO,jO,  wvins 
0.00  40  -80  0.2 

7)  option  for  snapshot  display:  component,AGC,begin  6  end  column 
"m"  1 
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1  500  150012X2 

XI  00  1  2001  2  IX 

XI  00  150012X0 

XI  00  X120X22X 
XI  00  XI  1X1  2  IX 
XI  00  1X1  1X1  2  IX 
1101  X  1X1  1X1  2  IX 
1X1  X  1X1  2161  2  IX 
1X1  X  mi  2341 2  IX 
5)  4  of  JUftowoynp/n 
31 

0)  coeWfcwHi  of  aontoTK  (lirw  ImJCittno  2mZ’$,ntgXmalnin  MOOMor) 
0120X044  OX  OX 032  IX  IX  140  IX  IX  1X204  2X 
2XX22XX4  3X316  332  345  364  3X3X412 4X444  4X4X402 
OX  OX  OX  OX  OX  OX  OX  OX  OX  OX  OX  OX  OX  OX 
0X034  018002002002002002002002002002002002002002002 

10)  toX*  of  OBtaliong; 4  of  OBmUomlinop^tot 
XO  IX 

11)  marching  grid  t>Om^yoa)?  toX  grid  aridfh  (Mntt 
0500 

EOF 

4mmm  ntntd 
fd2<MFD>  5  FDjnsg 

4mmm  dOOlUlOlltOX  IMM  80/104 

UapkK  daUL0>  5  omr 
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The  kt2  genecales  the  following  error  message 

fD  codi  Mwston  SSMOSL  RSJ 
Compitd  m0t  gikU  SCO  600 
OiUdaakwd  500  200 

Hath  6  tmfionl  Mfitfoktot  m  .0100  .0100  .0010 

HotaoQanaoua  taoOti 
Input  lapogiwfilv  OMtTOP 
S£TTOP>lih-  66 

SETrOP>  OpgO  maa  (mmpM  mmy  50  aokmmt): 

00000  SI  64646464 
seTTOP>  idh(l,3i  (aamplad  auaiy  SO  eokatam): 

606666666638  2  2  2  2 
Boundaiy  Saga  0030 
OnitUaaok^aaln  SETTOP) 
gik^  land(ootunm  1)m  65 
gik^  jandfeotumn  kuym  1 
QRtD>  gaotogical  maM  taadjf. 

QMaampbdpar  SlbpSInodaa 
UmbdilWtf 


2: 

6 

6 

.0 

.0 

.0 

6  21.0  21.0  21.0  216 

23: 

.0 

6 

6 

.0 

.0 

.0  21.0  21.0  21.0  216 

44: 

6 

.0 

.0 

.0 

.0  21.0  21.0  21.0  216  216 

68: 

.0 

6 

6 

.0 

.0  216  21.0  21.0  216  21.0 

66:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21j0  21.0  21.0 

107:  21.0  21.0  21.0  21.0  21.0  21.0  21J0  21J0  21M  21.0 

126:  21.0  21.0  21.0  21.0  21.0  21.0  21J0  21.0  21J>  21.0 

140:  21.0  21.0  21.0  21.0  21M  21.0  21.0  216  21.0  21.0 

170:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  216  216  216 

101:  21.0  21.0  21.0  216  216  21.0  21.0  21.0  216  21.0 

MuSaU 


2: 

.0 

.0 

.0 

.0 

.0 

.0  21.0  21.0  216  216 

23: 

.0 

.0 

.0 

.0 

.0 

.0  21.0  21.0  21.0  21.0 

44: 

.0 

.0 

6 

.0 

.0  21.0  21.0  216  21.0  216 

85: 

.0 

.0 

.0 

.0 

.0  21.0  21.0  21.0  21.0  21.0 

86:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  216  21.0 
107:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  216  21.0  21.0 

128:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  216 

140:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0 

170:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  216  21.0  21.0 

191:  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0  21.0 

DansilyfUd 


2: 

.0 

.0 

.0 

.0 

.0 

.0 

25 

25 

25 

25 

23: 

.0 

.0 

.0 

.0 

.0 

.0 

25 

2.5 

25 

25 

44: 

.0 

.0 

.0 

.0 

.0 

2.5 

25 

25 

25 

25 

65: 

.0 

.0 

.0 

.0 

.0 

2.5 

25 

25 

25 

25 

86:  2.5  2.5  25  2.5  25  25  25  25  25  25 

107:  25  25  2.5  2.5  2.5  2.5  25  25  25  2.5 

126:  25  2.5  25  2.5  2.5  25  2.5  25  25  26 

149:  2.5  2.5  2.5  2.5  2.5  2.5  2.5  25  25  2.5 

170:  25  2.5  2.5  2.5  2.5  2.5  25  25  2.5  2.5 

191:  2.5  25  2.5  2.5  2.5  2.5  25  25  25  2.5 

Sguarad  P  vahdty  Held 
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2: 

.0 

.0 

.0 

.0 

.0 

.0  252  252  252  252 

23: 

.0 

.0 

.0 

.0 

.0 

.0  252  252  252  252 

44: 

.0 

.0 

.0 

.0 

.0  ^2  252  252  252  252 

85: 

.0 

.0 

.0 

.0 

.0  252  252  252  252  252 

ML-  25^  2S^  25.2  2S.2  2&2  2&2  25J  2&2  2SJ! 

107:  25.2  2SJ!  25.2  25.2  25.2  25.2  25.2  252  252 

125:  25.2  252  25.2  25.2  25.2  252  252  252  252  X.2 

140:  252  252  25.2  25.2  252  252  252  252  252  252 

170:  252  252  25.2  252  252  252  252  252  252  252 

101:  252  252  252  252  252  252  252  252  252  252 

Squtnd  S  vahei^  lUd 


2: 

.0 

.0 

.0 

.0 

.0 

.0 

54 

54 

54 

54 

23: 

.0 

.0 

.0 

.0 

.0 

.0 

54 

54 

54 

54 

44: 

.0 

.0 

.0 

.0 

.0 

8.4 

8.4 

54 

54 

54 

65: 

.0 
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.0 

.0 

8.4 
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54 

54 

54 
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8.4 

54 

54 

54 

54 

8.4 

54 

54 

8.4 

54 

107: 

54 

54 

8.4 

8.4 

8.4 

8.4 

54 

54 

54 

8.4 

128: 

54 

54 

8.4 

8.4 

8.4 

54 

54 

54 

54 

54 

140: 

54 

54 

8.4 

8.4 

8.4 

54 

54 

54 

54 

54 

170: 

54 

54 

8.4 

8.4 

8.4 

8.4 

54 

54 

54 

54 

101: 

54 

54 

8.4 

8.4 

8.4 

8.4 

54 

54 

54 

54 

(VslVp)“2 

2: 

.00 

.00 

.00 

.00 

.00 

.00 

.33 

.33 

.33 

23 

23: 

.00 

.00 

.00 

.00 

.00 

.00 

.33 

23 

23 

.33 

44: 

.00 

.00 

.00 

.00 

.00 

.33 

.33 

.33 

.33 

.33 

65: 

.00 

.00 

.00 

.00 

.00 

.33 

.33 

23 

23 

.33 

86: 

.33 

.33 

.33 

.33 

.33 

.33 

.33 

23 

23 

.33 

107: 

.33 

.33 

.33 

.33 

.33 

.33 

23 

23 

23 

.33 

128: 

.33 

.33 

.33 

.33 

.33 

.33 

23 

23 

.33 

.38 

140: 

.33 

.33 

.33 

.33 

.33 

.33 

23 

23 

.33 

.33 

170: 

.33 

.33 

.33 

.33 

.33 

.33 

23 

23 

23 

.33 

191: 

.33 

.33 

.33 

.33 

.33 

.33 

.33 

.33 

.33 

.33 

Selactad  wa¥e  type  m  7 
0.40-90  .200000003 
m  component  wil  be  displayed 
AGCOag^  1 

Display  leK  right  (column),  start  snd  (dl),  lop,  boaom  (roar) 
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1  5001 

2  202 

301 
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2000 

2  202 

701 
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901 

1261 
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901 
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1201 

1561 

2  100 

1101 

0 

1381 

1801 

2  100 

1301 

0 

1681 

2161 

2  100 

1501 

0 

1941 

2341 

2  100 

31  seismometers 

12  28  44  60  76  92  108  124  140  156 
172  188  204  220  236  252  268  284  300  316 
332  348  364  3W  396  412  428  444  450  476 
492 

66666666666666666666 
666666665034  18  2  22 
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222222g222 

2 

AWiMraltoM.  100  sttpa  ptr  bnak 
Vacfortumofdlafilaeanmnia 
Oispk^llagm  6 

Plane  P  wma  (Ohnaka  or  Oeuaaian  ahapa) 
body>  P  vahelly  near  grid  boOomm  5.02 
bo(fy>  wavalat  apan  (km)-  .318 

bo(H>  Mlwaaalnnt atarta muaMy al  79mw 
bod^  wanmlal  haa  z-lang8i  nugh^  31  pokaa 
Mdent  Angle:  .OOdag 
Boundary  Flaga  1 120 

(lape  10)  dataO  opanadi 
FARPSInlaiiaoeatdapdtmldt/2~  101 

Snapahol  S1000:  now,loop.kauJauJdaf^  111  SCO  2  202 
Snapahol  S1001:  noar.k)op,kaujau,idepam  101  1  1  SOO  2  202 

S^ttnograma  output,  loop  101 

Snapahot  S1002:  neat,loopJtiuJauJdepam  201  1  1  500  2  202 

Sakmogtama  output,  loop  201 

Snapahol  S1003:  noar,loopJmijau,ldapam  301  1  1  SOO  2  202 

Sakmo^^tMa  output,  loop  301 

Snapahot  S1004:  noar,loopJaaiJau,ldepam  401  1  1  500  2  202 

Saiamogratna  output,  loop  401 

Snapahol  S100S:  naaf,loopJtaujauJdapam  SOI  1  1  800  2  202 

Salamogtama  output,  loop  801 

Snapahot  S1006:  na»tJoopJaujau,idepam  801  1  1  SOO  2  202 

Saiamograma  output,  loop  801 

Snapahot  S1007:  noar,loopj(aujau,ldepam  701  1  1  SOO  2  202 

Seiamogtatna  output  loop  701 
Saiamogram  and  anapahot  /Im  doaad 
WavaSalda  backup  dona  at  hataalon  701 
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wMiiatfjrw  “movm“ 


Sttp  3:  QMMratIng  Snapshots 

Bmmm  Mmfih  toift  to  plat  out  ana  'otato 
talnonomaloh 

$—  pnpan  kba!  for  bgurat 
eat«BOF>Ubal 
ISikgnmp,  Ckmulmt  pkrmr  wnm 
OM $00*202.  tbcmlOrn.  dklma 
HaKapmoa:  Vpm&fOkrnla,  V*m2MOMa 
EOF 

dmiga  Ota  labal 

sat  Ub^  S1000  S1001  S1002  $1008  ) 
fbraaehl(1 234) 
ax$klsffq«  EOF»  Junk 
4.6d 
3 

rLabaf 
1 

wq 
EOF 
end 

a— pid  out  snapshots  with  “moida" 
movia  -ar  J  3  -m  4  $kk>  3  Error 

Figure  10  shows  the  propagation  of  a  normaly  inckJenl  planar  Gaussian  P  wave  through  a  modei 
with  a  45*  ramp  on  the  top  of  the  grid  and  von  Neumann  (/.e.  O-slope  or  symmetric)  boundary  conditions 
on  both  sides.  The  homogeneous  mecflum  has  the  compressiortal  and  shear  velocities  of  5.02  and  2S98 
krrys,  respectively.  Shadkig  in  each  snapshot  is  proportionai  to  the  displacement  amplitude  -flFTW. 

The  accuracy  of  the  results  was  tested  empiricaly  by  comparing  results  obtained  from  two 
separate  fimteKfifference  nins.  In  the  first  nin.  wo  adopted  grid  spadngs  Ax  =  Az  -  0.01km.  In  the 
second  mn,  we  used  one-half  the  vertical  spacing  but  dO(4)led  the  vertical  exaggeration  at  discretization. 
Thus  two  algorfthms,  B  and  A  in  Section  2.3,  were  used  separately  in  these  two  finite-difference  runs. 
Excellent  agreement  was  obtained,  iiKlcating  that  the  two  algorithms  A  and  B  are  consistent  and  that 
both  solutions  are  good  approximations  to  the  true  solutions,  by  Lax's  equivalence  theorem  (Smfth, 
1978;  MHcheli  and  Griffiths,  1980).  The  apprtqxiate  P-SV  conversions  and  the  reflections  and 
diffractions  satisfying  SneO's  law  artd  Huygen's  principle  are  dearly  visible  in  the  successive  snapshots 
taken  every  0.1  second,  as  shown  in  Figure  10.  It  is  easy  to  verify  the  first-order  accuracy  in  spatial 
increment  of  our  one-sided  explicit  representation  of  the  free-surface  boundary  conditions  in  this  case. 
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acMT(v~9Mr*«  Moaec.  i  >oe 

XTIeltafll-  IlMI;  ZTIellWk-  H 


FDalap  too 

8aRT(v*‘a*ir*2)  .ioosec,  i  boo  i  202 

XTIcHMkB  tUnX;  ZTIeilaifcB  1  UnM 
■OOOOOfctOO  .20002E«01 


FDtl^  200 

SCMIT(V**2«H**2)  J006EC,  1  BOO  1  202 

XTteMario  1  UnH;  ZTIeltario  1  UnK 
JKX>0B6*00  .12B03e4O1 


FOatop  200 

8QRr(V**24H**2)  J008EC,  1  BOO  1  202 

XTIcMaili-  1  UnK;  Z  Tie  Mario  1  UnR 
MOOOE4OO  .I86ISE4OI 


45-deg  ramp,  Guassian  planar  wava 
Grid  500*202,  dxsdzsIOm,  dtsima 
Half  Space:  Vpr5.02km/s,  Vss2J98l(ni/a 


Figure  10 


step  4:  Plotting  Synthetics 

Figure  1 1  shows  the  vertical  components  recorded  at  31  grid  points  on  the  free  surface  for  630 
marching  steps  in  time,  with  4t  «  0.001  sec.  The  sensors  are  0.16km  apart  in  X  directioa  The  numeri¬ 
cal  label  (at  the  right)  irKiicates  the  distance  from  the  corresponding  sensor  to  the  nearer  comer  point. 

Reflections  from  the  edge  would  begin  to  contaminate  the  wavefield  after  0.6  sec,  due  to  the  symmetric 
boundary  conditions.  This  figure  is  generated  enth  a  utiitty  program  “arrayplt"  using  the  foRowing  script. 

Note  that  we  have  two  more  utility  programs  “trp»90“  and  “trpit"  which  read  exactly  the  same  input 
parameter  file  to  generate  different  types  of  plots  (e.g..  Figures  14  and  16). 

* 

set  nononuHch 

setenv  PSLANDSCAPE 

*—  rotate  the  pht  by  90  degrees 

set  DBm(  VERT  HOfV  ) 

(oreacb  1(1) 
cat«  EOF>  input 

1)  job  type  (oddm^.w,  evetM>ASCU,  negabvemmokaep  error  messages) 

2 

2)  database  name  (char*40): 

$DBl$ij 

3)  sampling  rate 
1000.0 

4)  a  o/  traces  to  be  ptoaed 
31 

5)  give  seismogram  Indices: 

1  23  4  56  78  9  10  11  12  13  14  15  16  1718  19  2021  2223  24  2526  2728  29  30  31 

6)  give  staring  point  of  each  trace: 

111111111111111111111111111111111111 

7)  give  signal  window  of  each  trace: 

630630630630630630630630630630630630630630630630630630630630630630630630630630630630630630 

8)  turn  on  de-spiker  (>0^>^ye8)?  turn  on  9fP/o-4e  detector  t>0^yes) 

00 

9)  AGC(0)  or  normalized  by  lmax(1),  constant(2),  specific  lrace(3)) 

2 

9-2)  give  norrtaUzing  constant 
0.30 

10)  black  out  negative  part?  (>0=^yes)  tic  maria?  (>0  yes) 

10 

1 1)  reduced  time  scale?  (>0  yes)  group  velocity  (dummy  if  .ie.O) 

0  10  0 

12)  decimation  option 
1 

13)  title 

45-deg  ramp,  grid  SOO'400,  dxmlOm,  dt=.001s,  plartar  P  wave 

14)  give  channel  ID's  or  distances  (as  many  lines  asS  of  traces) 

-2.08km 

1.92km 

-1.76km 

-1.60km 

-1.44km 
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•IJBkm 

-U2km 

•0.961m 

•O-BOm 

-0.641m 

•0.461m 

-0.32km 

■0.161m 

eomtrl 


ramp 

tamp 

oomar2 

0.16km 

0.321m 

0.46km 

0.64km 

0.60km 

0.96km 

1.12km 

IJOkm 

1.44km 

1.60km 

1.761m 

1.92km 

2.06km 

EOF 

anaypk  <  input  >  Satror 
and 


it  is  interesting  to  note  that  for  an  incident  P  wave,  the  bottom  comer  of  the  ramp  (or  step)  is  an 
efficient  scatterer  of  P  waves,  whereas  the  top  comer  of  the  inclined  surface  (at  45*^  produces  a  strong 
scattered  Rayleigh  field  {cf.  Figures  10  and  11).  The  bottom  comer  has  the  smallest  peak  ampHtude 
(cf .  the  rightmost  column  of  Figure  11),  whereas  ttie  top  oomer  has  the  largest  peak  ampHtude. 
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CO'^IO<DI^OOO>Ot- 

CMCMCMCJCVICMCMCOCO 


45-deg  ramp,  grid  500*400,  dx=dz=10m,  dt*.001s,  planar  P  wave 


6.  EXAMPLE  2:  RAYLEIGH  WAVE  SCATTERING  AT  SEOMENTARY  BASIN 


The  relative  importance  of  near-surface  effects  versus  deep-crustal  effects  on  relative  ampiiludes 
can  be  examined  by  modeino  in  laterally  heterogeneous  models.  Numerical  experiments  are  particulariy 
suitable  to  address  the  effects  of  variable  cnjstal  thickness,  njgged  mountainous  reHef,  sedknert  thick¬ 
ness,  and  other  stnjctural  heterogeneities  on  Lg  and  Rg  ampBudes.  Such  studies  would  oompiemeni 
the  empirical  studies  by  providing  better  interpretations  and  better  insigN  of  the  underlying  physics.  As 
an  illustrating  experiment,  here  we  give  an  example  to  demonstrate  how  LFD  calculations  can  be  utilized 
in  explaining  and  quantifying  the  Rg  /Lg  blockage  as  wel  as  the  higher  Sq/S„  coda  level  due  to  strong 
Lg  -to-Sn  and  Rg  -to-Sg  conversion  at  the  pinched  attenuating  sediment  layers. 

The  initial  waveform  is  chosen  so  that  the  vertical  displacement  W  on  the  flat  tree-surface  is  a 
Ricker  wavelet.  This  wavelet  has  been  adopted  frequently  in  finite-difference  simulations  becssjse  K  is 
well  localized  in  both  spatial  and  wavenumber  domains  (Boore.  1972;  Munasinghe  and  FameM,  1973; 
Fuyuki  and  Matsumoto.  1980;  Fuyuki  and  Nakano,  1984;  Levander.  1985).  To  avoid  grid  dtepersion,  we 
chose  1.2135  km  -  48.5  Ax  as  the  dominant  wavelength  of  the  incideit  Rayleigh  wave  packet  in  a 
homogeneous  portion  of  the  medium  with  a  P-wave  velocity  4.57  km/sec  and  Poisson’s  ratio  0.35. 
Absorbing  boundary  conditions  (Clayton  and  EngquisL  1977;  Emerman  and  Stephen,  1983)  are  adopted 
on  the  sides  and  bottom  to  suppress  the  artificiai  reflections  from  the  sides  of  the  grid.  Note  that  the 
quasi-transparent  boundary  conditions  do  allow  the  wave  to  disappear  into  the  sides  and  bottom  of  the 
grid. 

Step  1:  2-0  Model  Generation 

The  utility  program  “Iymed3"  generates  a  2-dimensional  medium,  of  which  the  material  properties 
are  constant  within  subregions  separated  by  (piecewise  linear)  polygons.  It  reads  the  coordinates  of  the 
vertices  of  polygons  and  fills  the  medium  above  the  polygon  with  prespedfied  material  properties.  By 
repeating  this  procedure  several  times,  this  routine  generates  a  2-dimensional  velocity  rrwdel  for  the  FO 
simulations.  “tymedSa"  and  “fymedSb"  produce  the  corresponding  graphic  presentation  of  the  medium 
(as  a  PostScript  output  file).  The  only  difference  between  these  two  codes  is  that  "tymedSa'Tit^  in  each 
layer  with  shaded  pattern,  whereas  “lymedSb"  fills  in  each  layer  with  a  predefined  pattern.  The 
PostScript  output  can  be  then  transmitted  to  laser  printers  directly. 

The  1st  line  of  the  input  file  specifies  the  label  (in  quotes),  which  is  followed  by  grid  size,  grid  spac- 
ings,  the  velocity  parameters  of  the  half  space,  as  weH  as  the  number  of  the  polygons.  ‘lymecO"  and 
“lymedSa"  then  read  exactly  as  many  polygons  as  defined.  Each  polygon  begins  with  the  number  of  ver¬ 
tices  in  this  particular  ooiygon  as  well  as  the  material  parameters  for  the  medium  above  this  polygon. 
The  material  parameters  are  typically  the  compressional  velocity,  shear  velocity,  and  the  density,  but 
they  could  also  be  the  Lame  constants.  These  parameters  are  not  required  for  the  graphic  routine 
"lymedSa". 

• 

set  nonomatch 

«— job  1:  njn  tymedS'  to  generate  medium  models  (Vp,  Vs,  fVto) 

«—  job  2,3:  run  1ymed3alb'  to  generate  PostScript  plot  of  the  model 
set  job=(  lyntedS  Iymed3a  lymedSb) 
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hmtehi(S) 

<M«l>kK30M 

800  500  (glmmadUa^kmJO^ 

0.005  0025  (gim  sptOal  ipmliio:  tt  kmf 

8.150520  (glmVp^V8J0iolnht2$pm») 

8  fhof/f  in8trtuP9tUtgtf9f 

24072.84  2.5  (t  cf  fioim at hOulmet  1:  Y^^ijtho atom  10 

1  122  (latarnmaHmihem  1) 

800  122  (2Hd  aamn  al  Indrtaea  1} 

45.001.80225  (t  at  palm  al  kmrtaoB  2  VP.^Rho  abaaa  10 

180  2  (lot  aaOn  at  Imriaea  2 

200  02  (2nd  vanmt  at  kttartaaa  2! 

520  82  (....) 

440  2  (...) 

4  2.140  1.1452000 
180  2 
220  52 
580  52 
440  2 

4  10400.840  1000 
100  2 
240  12 
580  12 
440  2 
/ 

^ol)t$0<UOD.IN 

and 


Program  ‘lymodSb'’ gononlm  the  foKwving  output  meeeage  as  wel  ae  a  PoetSortpr  formatted  fite 
which  can  be  directed  to  the  laser  printer  for  a  hard  copy  of  the  model  (Figure  12). 

LYUED5B.  950508  nsj 
"Sadbnantmy  Batin  (modal  B5)’ 
lymadStn  grid  atm:  OOP  800 
Iymad3t>  math  aridOit:  0250  .0250 

Iymad8>  htH-tpamt  pmtmtttn  (alpha,  halt,  iho):  0.150  5000  2000 

Iymad5b>  «  at  Mmlaett:  4 
imtritea  1  hot  Oaarliem 


vridi  Vp,Vls.fOio  abate  8dt  katritat: 
intttttea  2hat  dverOeaa 

4.5700 

28400 

2.5000 

wilh  yp,Vs,Rho  abate  84a  krittitet: 
immfaree  3hat  4vemom 

30000 

1.8000 

22500 

whh  Yb,Va,Rho  above  84a  kritrftea: 
katrftee  4  has  4  vetdott 

21400 

1.1430 

2.0000 

vridt  Vp.Vs,Rho  abate  d4t  kriariaot: 

1.3400 

.8400 

1.8000 

total  KZ(lm).  .15000£*02  .1^00E*02 
Xteala.ddx,ring8t~  .5353  .0083  5.0000 

Z$calB.ddz,2lnglh~  -0400  -.0080  -30000 

Inlarftce  1  (from  bottom) 

4  revised  Jordan  vertices: 
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vfi.vt.fnm  4Sfo  za40  zaoo 

1  0 

1  130 

600  120 

600  0 

Martaot  2(»nmbaoom) 

4  rmfiMd  Jordan  rafeaa: 
^Va.ftOm  XOOO  1.600  2200 
160  0 

200  60 

220  60 

440  0 

kmriaoa  2  (kom  boOom) 

4  rmiitad  Joidan  aortieaa: 
VfKVo,f»m  2.140  1.142  2jOOO 

160  0 

220  20 

260  20 

440  0 

mtmnm09  4  (wOm  OOm/fnJ 

4  faatoad  Jordan  aardoaa^ 
Vfi.Vt.Rhem  1.240  .640  1M0 

160  0 

240  10 

260  10 

440  0 
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'Sedimentafy  Basin  (model  E3)’ 

(only  uppormoat  layors  ara  ahown) 

0X=  .02S0Km ,  OZs  .(esOKin 
X  tic  mark  sZ  lie  nurks  1km 
Half  apaea  paramalara:  6.190  3J00  2.800 
imarfaca  1  from  boltom:  4.970  2JM0  2.900 
4  ravlaad  Jordan  vardoaa: 

(  1.  0)  (  1.  1201  (  600,  120) 
Inlarfaea  2fromboNom;  2.000  1M0  2.290 
4  ravlaad  Jordan  varUcaa: 

(  160  ,  0)  (  280  ,  60)  (  320  ,  60) 

imarfaca  3  from  bottom:  2.140  1.143  2.000 
4  ravlaad  Jordan  vatdcoa: 

(  160,  0)  (  220,  30)  (  380,  30) 

imarfaca  4  from  bottom:  1.340  .640  1.8(M 

4  ravlaad  Jordan  vartieaa: 

(  160,  0)  (  240,  10)  (  360,  10) 


600,  0) 

(  440,  0) 

(  440,  0) 

(  440,  0) 


Figure  12 


SMp  2:  Running  FD 

aat  nonomtteh 
ett«eOF>  kput 
grill  dhmmion  kmrjih 
900900 

KXtpaelHgtol9mgririm&ah9lmipoal»paeing:dKdc^ 

0.0280.0950.0029 

IwmogamoiM  or  htlangmmouM  (01.21 
2 

InVp 

kiV* 

wwwtO 

topography  modoi  kiTOP 
norm 

ChOOM^  IMMW 

4 

alpha.  OJO.  wrina 
990.009020  1.2198 

9»  op9on  tor  anapahot  rOapky:  oompormrajlQGJbaghi  9  amt  ookam 

"hr  1 

1  900  1  4001 22  200 
1  0  1  9001  22  900 

901  0  991  721  1  100 

401  0  901  991  1  100 

401  0  001  1291  1  100 
401  0  1201  1991  1  100 

101  0  1991  1901  1  100 
901  0  1991  2191  1  100 

401  0  1041  2941  1  100 

ataOona  for  raoonOng  aahmograma 
90 

eooninalaa  cf  aanaora:  (2na  1  m  JOa.  tna  2  m  Za) 

040090  120  190  200240290  920990400  440490  820090000 
040090  120  190  200240200820990400440400820990900 
002002002002002002002002002002002002  002002002 
200200200200200200200200200200200200200200200 
loM  number  of  hararioru;  OaraHonatanapahot 
4001200 

mardOng  grid,  total  aba 

0900 

EOF 

fd2<  input>  9  FD2jnag 
fdaph  <  data.0  >  9  Error jnag 
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step  3:  Generating  Snapshots 


Lg  and  Rg  are  susceptK)le  to  sedknenlary  basin  stnjctures  and  other  types  of  lateral  heterogeneMy 
such  as  rnoumainous  topography.  The  scattering  at  the  edge  of  a  pinched  attenuating  bas^  could  be 
very  severe,  as  illustrated  in  these  snapshots  of  the  displacement  wavefield.  Snapshot  3  irxScatM  that 
the  incident  Rg  generates  higher  modes  in  the  baan,  folloMred  by  slower  fundamertal  mode  Rayleigh 
waves,  as  well  as  the  converted  body  waves.  Note  that  a  large  buk  of  the  converted  S  wave  energy 
tunnels  underneath  the  basbi.  In  snapshot  5,  ancHher  strong  Rayleigh-to-S  conversion  at  the  edge  of  the 
basin  is  dearly  vistt)le.  This  LFD  example  can  be  scaled  to  explain  the  Lg  blockage  by  Barents  Shelf  as 
well  as  the  high  S„  coda  level  observed  at  NORESS  due  to  strong  Lg  -to-Sn  conversion  reported  by 
Baumgardt  (1991). 

* 

set  nonomalch 
cat«  />  me 
4 
2 

1  122 
600  122 

4 

160  2 
260  62 
320  62 
440  2 
4 

160  2 
220  32 
380  32 
440  2 
4 

160  2 
240  12 
360  12 
440  2 

I 

cet«  />  Utbel 
Sedimentary  Basin  (model  E3) 

600‘S00  [dxm^.dt~2.Sms] 
fig:  1.213Skm,  2Hz 
I 


S—  0,1,3,6.10  sec 

setkk  =  (  S1004  S1014  S1034  S1064  S1 104  ) 
foreach  1(12345  ) 

«—  change  label 
exSkkffi]  «  EOF»  junk 
4,6d 
3 

r  Label 
1 

wl 
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9 

EOF 

ana 


»— plat  out  anapaMa 
moMto  •iS-HJHa-mS$iat>tanor 
§movia-t-i3-tuaa-mS$kk>aanor 
ameuia -v -i  2 -I  LSa -m  8  $U(>  3  airor 
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FDalap  0 

aairT(v~2+H~3)  omosec,  i  aw  i  aoe 
XTkltarti>  lUna;  ZTIcMailis  1UM 
jioooee4ao  .i2vrK4«i 


FOalw  MO 

SQilT(V**2«H**a)  1.0006CC,  1  tW  1  3M 

XTteltoika  IlMI;  zncllatks  1  UnM 
W000E«00  .77114E«M 


FDalw  laOO 

aQivT(v**s«H~3)  ajMoaec,  i  aw  i  aw 

XTIellMkB  1  UnM;  ZneUafkB  lUhX 
jMOOoe+oe  jiis3E*oo 


roaiw  Z400 

SQI(T(V**2«H**2)  aJMOSCC,  1  600  1  200 
XHeltailca  lUnX;  ZTteHailo  lUhM 
MOOOC4OO  ja979C*00 


FOalap  4000 

SQIfT(V~2«H**2)  10M08EC,  1  600  1  200 

XTIeHartis  1  UnX;  ZTIcMaifcs  lUnK 
•OOOOOE+OO  .6S667E-01 


Ssdimentary  Basin  (model  E3) 
600*500  [dx=25m,dts2.Sms] 
Rg;  1.2135km,  2Hz 


Figure  13 


step  4:  Plotting  Synttiellcs 


tcript  for  '«no«plt'  or  ‘kpMO':  to  plot  out  m  bundt  otwoeoo 
set  nonomtPch 
satanv  PSLANt^CAPE 
set  De~(  VERT) 
foreachif  1 ) 


otS«  EOF>  Input 

1)  job  type  (oddmu>.w,  euenmu>ASCU.  negetyum^keep  amr  messegee) 

2 

2)  detabase  name  lchai‘40): 

soem 

3)  aemptng  rale 
400.0 

4)  a  of  tnoes  to  be  plotted 
IS 

5)  give  seismogram  Mhas: 

1  234S67891011  1213  14  IS 

6)  give  starting  point  of  each  trace: 

111111111111111111111111111111111111 

7)  give  signPI  wMoer  of  each  trace: 

SOOO  SOOO  SOOO  SOOO  SOOO  SOOO  SOOO  SOOO  SOOO  »J00  SOOO  5000  sooo  sooo  sooo  sooo  sooo  sooo  sooo  sooo 

8)  turn  on  de-sptter  (>0~^yea)?  turn  on  dtofoctor  ^O-eyes) 

00 

8)  AGC(0)  or  normetaed  by  {maxfl),  oonstant(2).  spedBo  trace(3)] 

1 

9-0)  give  normalizing  constant 

1111111111111111111111111111111111 

10)  black  out  negative  part?  (>0^yes)  tie  marka?  (>0  «»  yea) 

00 

1 1)  reduced  time  scale?  (>0  yes)  group  velocity  (dummy  if  Je.O) 

0  10  0 

12)  decimation  ration 
1 

13)  title 

Rg;  sedknentary  basin  (model  E3) 

14)  give  channel  ID's  or  dislanoes  (as  marry  Ones  as  •  oftaces) 

•3km 

-2km 

-1km 

0km_B 

1km~B 

2kn\_B 

3km_B 

4kmJ3 

SkmJB 

6km_B 

7km~B 

+  1km 

*2km 

+3km 

*4km 
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EOF 

tpUO<lnpul>  *mnr 
m  hamper 
and 


7.  EXAMPLE  3:  RIPPLE  FIRING 


A  major  issue  for  the  Non-Prolieration  Treaty  is  the  discrimination  large  chemical  explosions 
from  possible  clandestine  or  small  nuclear  tests.  Unless  discrimination  is  posst)le,  the  rumerous  mMng 
blasts  could  provide  ample  opportunity  for  oonceaing  clandestine  tests.  Rippled-fired  explosions  are 
commonly  used  to  figment  rocks  dursig  quarry  arKf  open-pit  mining.  The  periodicity  inherent  in  the  rip¬ 
ple  firing  could  produce  a  seismic  reinforcement  at  the  frequency  of  the  delay  between  shots  or  rows.  It 
has  been  suggested  that  the  convolution  of  a  sirxjle  explosion  wkh  a  comb  function  of  variable  spacing 
and  variable  amplitude  can  be  used  to  model  the  distinctive  signature  of  ripple  firirx)  (Stunp,  1988; 
Anderson  and  Stump,  1989;  Smith,  1989;  Hedfin  et  ai,  1990;  Chapman  ef  al.,  1992;  Reamer  at  a!.. 
1992;  and  many  others).  BaumgardI  and  Ziegler  (1988)  delicately  demonstrated  that  the  incoherent 
array-stack  spectra  can  be  used  to  identify  some  multiple  shots  recorded  at  NORSAR.  By  superposi¬ 
tioning  the  waveform  due  to  a  single  shot  with  proper  time  delay,  they  were  able  to  model  the  source 
multipiicty  under  the  assumption  that  the  spatial  spreading  of  the  shots  is  negligble  with  respect  to  the 
distance  to  the  receiver.  The  work  by  Stump  et  a/,  successfully  characterized  the  major  features  of  the 
wavefield  due  to  ripple  firings  at  near-source  ranges. 

There  are,  however,  some  other  wave  exckation  characteristics  of  ripple-fired  explosions  which  are 
not  predicted  by  such  spectral  or  waveform  superposition  approaches.  In  tNs  example,  the  linear  finite- 
difference  (LFD)  method  is  utilized  to  seek  some  ir^ights  into  the  ripple-fired  explosions. 

The  preliminary  results  suggest  that  ripple  firing  could  excite  directionally  dependent  Rg  and  S*. 
Thus  the  lack  of  Rg  may  not  be  always  indicative  of  a  deep  source.  Rather,  it  could  also  be  due  to  the 
shot  pattern.  However,  the  enhanced  Rg  in  the  forward  direction  of  ripple  firing  can  be  strongly 
attenuated  by  lateral  heterogeneity  and  surface  topography.  The  scattered  Rg  energy  could  then  couple 
into  the  crustal  waveguide  as  Lg  and  other  phases  (McLaugNin  and  Jih,  1986, 1987;  Jih  and  McLaugh¬ 
lin,  1988).  Since  such  scattering  mechanisms  are  commonly  present  in  many  quarry  sites  or  mines,  it  is 
not  surprising  that  the  directional  enhancement  of  Rg  may  not  be  always  obsenrable.  The  spaU  could 
obscure  the  azimuthal  dependency  of  Rg  as  well.  Previous  LFD  modeling  studies  (McLaughlin  and  Jih. 
1986,  1987)  suggest  that  the  near-source  Rg  -to-S  scattering  is  usually  stronger  than  that  of  Rg  -to-P, 
which  could  provide  a  plausible  explanation  of  why  quarry  blasts  and  mining  blasts  should  discriminale 
less  well  from  earthquakes  than  would  contained  nuclear  explosions.  Further  quantitative  analysis  along 
this  line  could  be  very  useful. 

Step  1:  Model  Generation 

Suppose  we  want  to  test  7  ripple-fired  explosions  in  a  model  with  a  trapezoidal  topography.  The 
whole  grid  has  a  dimension  of  800  by  600  nodes  with  the  topography  embedded  in  the  central  portion. 
The  following  sample  script  generates  the  topogrsqphy  file  for  the  nominal  model  in  which  we  set  off  the 
ripple  fires  from  left  to  right.  The  topography  starts  with  a  flat  segment  of  350  grids,  a  45**  ramp  of  30 
grids  high,  a  plateau  of  100  grids  long,  a  -45**  ramp  of  30  grids  long,  and  followed  by  a  flat  segment  of 
290  grids.  To  make  sure  that  the  topography  is  i;hown  in  the  central  portion  of  the  snapshots,  we  will 
extract  the  wavefield  from  columns  231  to  630.  The  topography  for  the  flipped  model  would  have  a 
reversed  order,  and  we  will  extract  from  columns  171  to  570. 
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• 

car«  />  kiTOPjigM 

32 

000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
00000000  00  00  0000  00  0000000000000000  000000  00  000000  0000000000  00 
00  00  00  00  00  00  00  00  00  00  00  00  00  00  00  OO  00  00  00  00  00  00  00  00  00  00  00  00  00  00 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
0000000000000000000000000000000000000000 

01  0203  04  05  06  0708  09  10  11  12  13  14  1518  17  18  10  20  21  22  23  24  25  26  27  28  20  30 
303030303030303030303030303030303030303030303030303030303030 
303030303030303030303030303030303030303030303030303030303030 
303030303030303030303030309090303030303030303030309030303030 
90903090909090909030 

9020  28  2726  2524  23222120  10  18  17181514  131211  1000  080706  0504  080201 
000000000000000000000000000000000000000000000000000000000000 
00  00  0000  000000  00  000000000000000000  000000  00  00  00  00000000  00  00  00 
000000000000000000000000000000000000000000000000000000000000 
00  00  0000  00  0000  00  00  00  000000000000  00  000000  00  00  0000  00  0000  00  00  00 
000000000000000000000000000000000000000000000000000000000000 
0000  00  00  00  000000000000  000000000000  060000  0000  0000  000000000000 
ooooooooooooooooooaooooooooooooooooooooooooooooooooooooooooo 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 
oooooooooorooooooooooooooooooooooooooooooooooooooooooooooooo 

0000000000  i.'^JOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 
000000000000000000000000000000000000000000000000000000000000 
000000000000000000000000000000000000000000000000000000000000 

I 
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step  2:  Running  FO 

For  the  nominal  model,  the  7  eouioes  are  located  at  <385, 5>,  <:390,S>,  <395,5>,  <400,5>,  <405,5>, 
<410,5>,  and  <415, 5>.  Except  the  first  shot,  which  is  detonated  at  tbne  step  0,  the  remaining  6  explo¬ 
sions  are  detonated  at  20  time  steps  apart  (cf.  the  parameter  file  b^ow).  For  the  flipped  model,  we 
reverse  the  order  to  set  off  the  first  shot  at  time  step  0  at  grid  point  <415.5>.  In  the  nominal  model,  we 
show  the  wavefield  between  columns  231  and  630,  whereas  for  the  flipped  model,  we  show  the  grid 
between  columns  171  and  570.  Thus  once  we  flip  the  model  back,  both  mns  will  be  showing  physically 
the  same  portion  of  the  wavefield. 

tor  nominal  modal 
sat  nonomatch 
»naad:  DRIVERS 

*  ganerata  Input  panmatar  Raforkl 
cat«  EOF>  InFD 

1)  grid  dmansion  kw,kh 
600600 

2)  x,z  spadngs  of  6>a  grid  mesh  6  tamponl  spacing:  dx,dt,dt 

0.01  0.01  0.001 

3)  homoganaous  (0),  hetamganeous  (1^ 

0 

InVp 

InVs 

InRho 

4f  topography  Oe 
inJOP 

5)  choose  amve  type 
9 

S-1)  drner  We 
DRIVERS 

S-2)  how  many  ripples 
7 

020390S 
040  3955 
050  4005 
050  4055 
100  410  S 
120  41SS 

6)  alfrita.  iOjO,  wrins 
20.00  36SS  0.2 

7)  option  for  snapshot  display:  componant^GC.bagin  6  end  column 
“V-  1 

1  600  1S20122  240 

231  630  1  2201  22  120 

201  00  1S0012  600 

301  00  601  200022S0 
701  00  901  1261  2  100 
901  00  1201  1S61  2  100 
1101  00  1361  16012  100 
1301  00  1661  2161  2  100 
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1801  00  194128412100 

8) 4«fm4$mo9mplM 
10 

9)  eootdkmrnt  of  atnton:  (Kh0  ImX'iJkm  mtteri 

001  100200900400800800  700000001  100  200  900 400  800  800  700  800 
092  092  092  092  OV 092  not  002  OOKKIt  m  OOP  OOt  002  002 

1(9  toml9  dOtnOom:*  oUtmatonOmaoahet 
4201  100 

11)  nmehlng  grid  (>Om^ymt)f  tdO  giU  niUOt  Oaalmd: 

4001200 

EOF 

U2<lnR)>tFDjfuo 
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Step  3:  Generating  Snapstiota 

In  the  input  parameter  file  for  “ki2“.  we  spedfied  two  subgrids  for  snapshot  output  {cf.  the  dtecitf- 
sion  in  Step  2).  The  foHowkig  is  a  sample  script  for  generating  the  snapshots  of  the  central  portion  only. 
Note  that  the  topography  radiates  the  fundamental  RayieicMrave  energy  into  body  waves  of  much 
lower  frequency,  as  shown  at  0.5  and  0.7  sec.  Part  of  the  Rayleigivwave  energy  is  trapped  between 
the  two  ramps  (cf.  0.6  and  0.7  sec)  and  eventualy  scattered  into  the  coda  waves.  The  topogrt^thy 
attenuates  the  Rayleigh  wave  in  both  directions  dramaticaly  and  produces  a  more  complex  waveform. 

t 

sat  nonomatch 

$—  2  t^pes:  Right  (hr  nominal  model),  Lett  (hr  thpad  modal) 

sat  typam(  Right ) 


ch«EOF>Lat>al 
7  shatom  ripph-firad  Shota 
TrapaiokU  hpography 
y/p~S.02limls.  Vs~2.80ekmts 
EOF 

eat«  />  Imp 

3855 

3905 

3955 

4005 

4055 

410  5 

4155 

I 

ll($typa  » telf  )hen 
satmovam(  170 ) 

cat«  l>  mpa/f.l 
program  Hipwf 

930304;  mp  over  character  Helds 
characler‘800  A,B 
characler^lO  label 
read(5,  '(a  10,2iS)')b4)elJong.mn 
wrile(6,  '(a  10.2S)')label,long/nn 
dohl.ll 
read(S,'(a)')A 
ijk’4nUnk(A) 
writ0(6,  '(aiyill:^) 
end  do 

100  raad(S.(a)'.end=200)A 

do  hi, long 
j=long+  1-i 
B(a)^(i:i) 
end  do 

write(6,  '(a)‘)B(1dong) 
go  ID  100 
200  stop 

end 
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MtoSmow 

$alklHSg001  8S009SZ006S2007S20M) 

hnmehi(12S4S) 
f _ e/MfUf  itbtl 

tlpmf<$U4$II>$K$ll 

•(M 

ep$um$KPI 

mtdK 

mm$9«BOF»Junk 

4.U 

S 

rUM 

1 

wq 

BOf 

and 


nnpa* 

cata<l>DO 

amk'ipriM-dmaaa.y<Tmp>  RUa 
I 

cihDO;miDO 

a— pht  out  anapahcta 

movia -t -i  S-r  (Via -mSfiy*  Error 
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VEHT  HOMKC.  ai  OO  1 
XTIelMiB  lUnX;  Z  He  Marti-  lUM 
.«Mi»iP0  j«0nc«oe 


PDaMp  MO 

vmr  t.Moioc,  mi  om  i  im 
XTIeMaifc-  1 1MI;  ZTk Marti-  1  UM 
-JM01t«00  A7MS*00 


PDala*  MO 

vmr  owoooocc;  zti  om  i  tM 
xneMaik-  IIMM  ZTIa Marti-  1  UnM 
■■MMMiM  ^OTOBtOO 


PDaMp  TOO 

vmr  O.TOOOCC.  mi  om  i  im 
XHeMa*-  lUnO;  ZTIa Marti-  1UM 
•JOOnEiW  .lOSTMMO 


POalap  000 

vmr  ojoooEC,  mi  om  i  im 

XHeMarii-  lIMt;  Z  Tie  Math-  lUnO 
•JSTWMtO  .lOIOOBiOO 


7  Shallow  rippla-flrod  shots 
Trapszoldal  topography 
Vps5.02km/s,  Vss2.898kin/s 


Figure  15 


step  4:  Merging  Synthetics 


As  specifled  in  the  sample  input  parameter  file  for  ‘id2’’  (page  53),  the  origkiai  grid  size  was 
defined  as  800x600.  As  soon  as  the  wavefront  gets  dose  to  the  edge,  the  program  ‘YioE?"  shifts  the  grid 
laterally  by  400  grid  points  {cf.  the  line  marked  “marching  grid”),  and  it  starts  to  save  the  synthetics  in 
another  buffer  called  “data.1”  (as  opposed  the  the  nominal  buffer  “data.0”).  Since  the  total  grid  width 
desired  was  1200  (-  800  +  400)  points,  only  one  shift  would  be  required.  At  the  end  of  the  LFO  run,  we 
need  to  merge  the  synthetics  aooorrflrig  to  their  physical  locations.  The  following  script  assumes  that 
the  buffers  “data.O”  and  “data.l"  are  stored  in  separate  directories  “Data.1"  and  “Data.2”,  respectively. 
The  demultiplexer  routkie  ‘idspH"  is  called  twice  to  decode  the  synthetics  in  these  two  directories.  Syn¬ 
thetics  associated  with  the  same  location  are  then  merged  (or  zero-flUed)  and  renamed. 

• 

set  nonomatoh 
tofeaeh  i  (  Lett ) 
fonach  i(  1  2) 
odSHDataj^ 
if  ( -a  data." )  then 

uncompraas  data.‘J^ 

Usplt  <  data*  >  A  Emr 
compiass  VEfTT 
rm  Etror  HOn* 

andd 

cdp/dtem 

end 

end 


eal«  i>  Z£ROj« 

0000000000000000000000000 

I 

cat  ZEROJS  ZEROJS  ZEROJS  ZEHO.ZS  >  ZEROJOO 

cat  ZEROJOO  ZEROJOO  >  ^R0_200  ~ 

cat  ZEROJOO  ZEROJOO  ZERO_26o  ZERq_200  >  ZEROJOO 

cat  ZEROJOO  ZEROJOO  >  ZEmjOO 

nnZ* 

a—mlgft 

@im1 

zcat  Lafl/DataJ/VERT.8ZLefltOataJVERT.4Z>  Z$i 
i  •${  +  1 

cp  ZEROJOO  ZSi ;  zcat  Lefi/DataJVERT.6.Z»  Z.$i 
@  i  =  $1  +  1 

cp  ZEROJOO  Z$i ;  zcat  LefllDataJVERT.eZ»  Z.$i 

#ss3  right 
@  f 

zcat  fegh»DataJ/VERT.8Z  Right Data_2/VERT.4.Z>  ZSi 
@  i  =  $i  +  1 

cp  ZEROJOO  ZSi ;  zcat  RghtDataJVERT.eZ  »  Z.Si 
@  i  =  Si  *  1 

cp  ZEROJOO  ZSi :  zcat  RightDataJVERT.8Z»  ZSi 
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mZERCy 


St«p  5:  Plolting  SynttMlics 

The  ooiTibined  synthetics  aie  then  ploled  out  with  yet  anottwr  plotting  routine  ‘trpff".  which  plots 
the  traces  in  two  columns.  In  this  case,  no  rotatino  of  the  plot  is  required.  Note  that  ‘trpt"  reads  the 
same  fonnal  of  input  parameter  file  as  do  “anaypT  ana  “ttpKSO". 

Brnrnrnieripthr'tpar  —  toplotoutmbimehcttaeu 

aatnononmleh 

MttDB^ZVERT) 

hnmOt  i(1 ) 


eat«  EOf>  kipul 

1)  job  OP*  (otkU^.w,  ooonmmoAaca.  mgatim  n  koof)  oaor  moougot) 

2 

2)  dMabOM  MHM  {idm^AOj: 

soem 

9)  oompang  mio 
1000.0 

4) »  et  trmooo  to  bo  lAMod 
10 

5)  gbm  toiomogrom  Mksoo: 

129994See« 

Ojglvmgboilngpointofoochtmeo: 

111111111111111111111111111111111111 

7)  gbm  tignol  toMoof  ofoaeh  bmoo: 

4000400040004000  4000400040004000  4000  4000 

8)  turn  on  OB-tptimr  (^Omm»ym)7  turn  an  SS)b-i>  dWBrtar  (>0¥>yaa) 

00 

8)  AQC(0)  ornommbadby  (imxffA  cwnfi(f2!l.  tpadOc  tnca(9)) 

2 

9-0)  give  nomwBzIng  oonattnt 

0.02S  0.02S  0.028  0.028  0.028  0.028 0.028 0j028  0.028  0.028  0.028  0028  0.0280.028  0.028  0.028  0.028 

10)  Uack  out  nogathn  part?  pO«^yaa)  lie  amke?  (>0  -•»  yea) 

00 

11)  radueed  tkna  aeale?  (>0  mmu  ya^  group  aalaeiljt  (dunrniy  8  Ja.O) 

0  100 

12)  decimation  option 
2 

19)titia 

Trapuzoid,  7  ahaOow  rippla-IM  blaata 

14)  giua  channel  ID'm  or  dislanaas  (aa  maiv  Knas  08  »  of  traces) 

-31m 

-Skm 

-7km 

-7km 

-7km 
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*Skm 

*-7lm 

*7lm 

*7lm 

EOF 
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